Es fascinante poder conocer ciertos aspectos del mundo del software libre y de código abierto orientado al análisis de datos espaciales, en está oportunidad quiero mostrarte un pequeño entrenamiento integrando dos softwares super potentes dentro del campo de la geocomputación, el primero es SAGA, el cual es un software de Sistema de Información Geográfica de código abierto y el otro es R, el cual es un lenguaje de programación libre e interpretado y orientado a objetos, este mismo presenta un gran ecosistema de paquetes para una diversidad de ejes temáticos, entre ellas podemos mencionar a la estadística, la estadística espacial, la teledetección, web mapping, machine learning, entre otros.
Para este entrenamiento tenemos que tener ciertos requerimientos como tener instalado SAGA y R, de preferencia se recomienda usar la IDE de Rstudio para tener una forma más amigable al momento de programar con R.
Se recomienda tener instalado la versión de SAGA 2.3 LTS hasta la 7.0 .
Si tienes instalado QGIS no es necesario instalar SAGA ya que viene por defecto dentro su ecosistema.
Por último, para integrar SAGA con R tenemos que instalar el paquete RSAGA (Alexander Brenning et.al,2018), posteriormente debemos identificar el directorio de trabajo donde está alojado la ruta de los módulos de SAGA.
libs <- c('RSAGA', 'tidyverse','raster',
'sp','sf', 'tmap', 'cptcity')
install.packages(libs, dependencies = TRUE)
library(tidyverse) # Pkgs para ciencia de datos
library(RSAGA) # Pkgs para integrar SAGA con R
library(raster) # Pkgs para manejar datos raster
library(raster) # Pkgs para manejar datos vectoriales
library(tmap) # Pkgs para elaborar mapas temáticos
library(cptcity) # Pkgs para obtener paleta de colores
Observación:
env <- rsaga.env(path = 'C:/Program Files/QGIS 3.10/apps/saga-ltr')
env <- rsaga.env()
## Search for SAGA command line program and modules...
## Done
env
## $workspace
## [1] "."
##
## $cmd
## [1] "saga_cmd"
##
## $path
## [1] "/usr/bin"
##
## $modules
## [1] "/usr/lib/x86_64-linux-gnu/saga"
##
## $version
## [1] "2.3.1"
##
## $cores
## [1] NA
##
## $parallel
## [1] FALSE
##
## $lib.prefix
## [1] "lib"
El proceso anterior nos ayuda a integrar SAGA con R, adicionalmente nos brinda una información acerca de la versión de SAGA con la cual estamos trabajando, la ruta de los módulos y alguna información adicional de interés.
libs <- rsaga.get.libraries(path = env$modules)
libs
## [1] "climate_tools" "contrib_perego"
## [3] "db_odbc" "db_pgsql"
## [5] "docs_html" "docs_pdf"
## [7] "garden_3d_viewer" "garden_fractals"
## [9] "garden_games" "garden_learn_to_program"
## [11] "garden_webservices" "grid_analysis"
## [13] "grid_calculus_bsl" "grid_calculus"
## [15] "grid_filter" "grid_gridding"
## [17] "grid_spline" "grid_tools"
## [19] "grid_visualisation" "imagery_classification"
## [21] "imagery_maxent" "imagery_opencv"
## [23] "imagery_photogrammetry" "imagery_segmentation"
## [25] "imagery_svm" "imagery_tools"
## [27] "imagery_vigra" "io_esri_e00"
## [29] "io_gdal" "io_gps"
## [31] "io_grid_image" "io_grid"
## [33] "io_shapes_dxf" "io_shapes"
## [35] "io_table" "io_virtual"
## [37] "pj_georeference" "pj_proj4"
## [39] "pointcloud_tools" "pointcloud_viewer"
## [41] "shapes_grid" "shapes_lines"
## [43] "shapes_points" "shapes_polygons"
## [45] "shapes_tools" "shapes_transect"
## [47] "sim_cellular_automata" "sim_ecosystems_hugget"
## [49] "sim_erosion" "sim_hydrology"
## [51] "sim_ihacres" "sim_qm_of_esp"
## [53] "sim_rivflow" "statistics_grid"
## [55] "statistics_kriging" "statistics_points"
## [57] "statistics_regression" "ta_channels"
## [59] "ta_compound" "ta_hydrology"
## [61] "ta_lighting" "ta_morphometry"
## [63] "ta_preprocessor" "ta_profiles"
## [65] "ta_slope_stability" "table_calculus"
## [67] "table_tools" "tin_tools"
## [69] "tin_viewer"
# Usamos la función lenght para indentificar el total de liberías:
length(libs)
## [1] 69
# Estos son los módulos para una libería específica ('ta_channels')
mod <- rsaga.get.modules(libs = 'ta_channels', env = env)
# Ahora vamos a crear una tabla con toda las liberías
# y sus respectivos módulos:
rsaga_mod <- rsaga.get.modules(libs = libs,env = env) %>%
bind_rows(.id = 'librerias') %>% as_tibble()
# Solo visualizamos las 10 primeras filas de toda la tabla
head(rsaga_mod,10)
## # A tibble: 10 x 4
## librerias code name interactive
## <chr> <dbl> <chr> <lgl>
## 1 climate_tools 0 Multi Level to Surface Interpolation FALSE
## 2 climate_tools 1 Multi Level to Points Interpolation FALSE
## 3 climate_tools 2 Earth's Orbital Parameters FALSE
## 4 climate_tools 3 Annual Course of Daily Insolation FALSE
## 5 climate_tools 4 Daily Insolation over Latitude FALSE
## 6 climate_tools 5 Monthly Global by Latitude FALSE
## 7 climate_tools 6 PET (after Hargreaves, Table) FALSE
## 8 climate_tools 7 Daily to Hourly PET FALSE
## 9 climate_tools 8 PET (after Hargreaves, Grid) FALSE
## 10 climate_tools 9 Sunrise and Sunset FALSE
dem <- raster('Insumos/DEM.tif')
plot(dem)
rsaga.fill.sinks(in.dem = 'Insumos/DEM.tif',
out.dem = 'Insumos/DEM_fill')
## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_preprocessor
## library : Preprocessing
## tool : Fill Sinks (Planchon/Darboux, 2001)
## author : Copyrights (c) 2003 by Volker Wichmann
## processors : 4 [4]
##
## loading: DEM
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## DEM: DEM
## Filled DEM: Filled DEM
## Minimum Slope [Degree]: 0.010000
##
## Save grid: Insumos/DEM_fill.sgrd...
# Visualización de nuestro DEM corregido
tmap_mode('view')
raster('Insumos/DEM_fill.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',alpha = 0.75,
palette = cpt(pal = 'ncl_topo_15lev'))
rsaga.slope.asp.curv(in.dem = 'Insumos/DEM_fill.sdat',
out.slope = 'Salida/pendiente',
out.aspect = 'Salida/aspecto',
out.cgene = 'Salida/curvatura',
unit.slope = 1,
unit.aspect = 1,
method = 'maxslope')
## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library : Morphometry
## tool : Slope, Aspect, Curvature
## author : O.Conrad (c) 2001
## processors : 4 [4]
##
## loading: DEM_fill
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Slope: Slope
## Aspect: Aspect
## General Curvature: General Curvature
## Profile Curvature: <not set>
## Plan Curvature: <not set>
## Flow Line Curvature: <not set>
## Method: maximum slope (Travis et al. 1975)
## Slope Units: degree
## Aspect Units: degree
##
## Save grid: Salida/pendiente.sgrd...
## Save grid: Salida/aspecto.sgrd...
## Save grid: Salida/curvatura.sgrd...
# Ahora procedemos a visualizar el mapa de pendiente
raster('Salida/pendiente.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'quantile',n = 7,palette = 'Spectral')
# Visualizando el mapa de aspectos
raster('Salida/aspecto.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',palette = cpt(pal = 'grass_gdd',rev = TRUE))
# Visualizando el mapa de curvatura general
raster('Salida/curvatura.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',palette = cpt(pal = 'grass_ramp'))
rsaga.get.usage(lib = 'ta_hydrology',module = 15)
## Search for SAGA command line program and modules...
## Done
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_hydrology
## library : Hydrology
## Usage: saga_cmd ta_hydrology 15 [-DEM <str>] [-WEIGHT <str>] [-AREA <str>] [-SLOPE <str>] [-AREA_MOD <str>] [-TWI <str>] [-SUCTION <double>] [-AREA_TYPE <str>] [-SLOPE_TYPE <str>] [-SLOPE_MIN <double>] [-SLOPE_OFF <double>] [-SLOPE_WEIGHT <double>]
## -DEM:<str> Elevation
## Grid (input)
## -WEIGHT:<str> Weights
## Grid (optional input)
## -AREA:<str> Catchment area
## Grid (output)
## -SLOPE:<str> Catchment slope
## Grid (output)
## -AREA_MOD:<str> Modified Catchment Area
## Grid (output)
## -TWI:<str> Topographic Wetness Index
## Grid (output)
## -SUCTION:<double> Suction
## Floating point
## Minimum: 0.000000
## Default: 10.000000
## -AREA_TYPE:<str> Type of Area
## Choice
## Available Choices:
## [0] absolute catchment area
## [1] square root of catchment area
## [2] specific catchment area
## Default: 1
## -SLOPE_TYPE:<str> Type of Slope
## Choice
## Available Choices:
## [0] local slope
## [1] catchment slope
## Default: 1
## -SLOPE_MIN:<double> Minimum Slope
## Floating point
## Minimum: 0.000000
## Default: 0.000000
## -SLOPE_OFF:<double> Offset Slope
## Floating point
## Minimum: 0.000000
## Default: 0.100000
## -SLOPE_WEIGHT:<double> Slope Weighting
## Floating point
## Minimum: 0.000000
rsaga.geoprocessor(lib = 'ta_hydrology',module = 15,
param = list(DEM = 'Insumos/DEM_fill.sdat',
TWI = 'Salida/Twi' ),
env = env)
##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_hydrology
## library : Hydrology
## tool : SAGA Wetness Index
## author : (c) 2001 by J.Boehner, O.Conrad
## processors : 4 [4]
##
## loading: DEM_fill
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Weights: <not set>
## Catchment area: Catchment area
## Catchment slope: Catchment slope
## Modified Catchment Area: Modified Catchment Area
## Topographic Wetness Index: Topographic Wetness Index
## Suction: 10.000000
## Type of Area: square root of catchment area
## Type of Slope: catchment slope
## Minimum Slope: 0.000000
## Offset Slope: 0.100000
## Slope Weighting: 1.000000
##
## Create index: DEM_fill
## catchment area and slope...
## pass 1 (39594 > 0)
## pass 2 (9653 > 0)
## pass 3 (4059 > 0)
## pass 4 (1920 > 0)
## pass 5 (950 > 0)
## pass 6 (484 > 0)
## pass 7 (246 > 0)
## pass 8 (130 > 0)
## pass 9 (64 > 0)
## pass 10 (32 > 0)
## pass 11 (18 > 0)
## pass 12 (5 > 0)
## pass 13 (2 > 0)
## pass 14 (0 > 0)
## post-processing...
## topographic wetness index...
## Save grid: Salida/Twi...
rsaga.get.usage(lib = 'ta_morphometry',module = 18)
## Search for SAGA command line program and modules...
## Done
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library : Morphometry
## Usage: saga_cmd ta_morphometry 18 [-DEM <str>] [-TPI <str>] [-STANDARD <str>] [-RADIUS_MIN <double>] [-RADIUS_MAX <double>] [-DW_WEIGHTING <str>] [-DW_IDW_POWER <double>] [-DW_IDW_OFFSET <str>] [-DW_BANDWIDTH <double>]
## -DEM:<str> Elevation
## Grid (input)
## -TPI:<str> Topographic Position Index
## Grid (output)
## -STANDARD:<str> Standardize
## Boolean
## Default: 0
## -RADIUS_MIN:<double> Radius
## Value range
## -RADIUS_MAX:<double> Radius
## Value range
## -DW_WEIGHTING:<str> Weighting Function
## Choice
## Available Choices:
## [0] no distance weighting
## [1] inverse distance to a power
## [2] exponential
## [3] gaussian weighting
## Default: 0
## -DW_IDW_POWER:<double> Inverse Distance Weighting Power
## Floating point
## Minimum: 0.000000
## Default: 1.000000
## -DW_IDW_OFFSET:<str> Inverse Distance Offset
## Boolean
## Default: 1
## -DW_BANDWIDTH:<double> Gaussian and Exponential Weighting Bandwidth
## Floating point
## Minimum: 0.000000
rsaga.geoprocessor(lib = 'ta_morphometry',module = 18,
param = list(DEM = 'Insumos/DEM_fill.sdat',
TPI = 'Salida/Tpi' ),
env = env)
##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library : Morphometry
## tool : Topographic Position Index (TPI)
## author : O.Conrad (c) 2011
## processors : 4 [4]
##
## loading: DEM_fill
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Topographic Position Index: Topographic Position Index
## Standardize: no
## Radius: 0.000000; 100.000000
## Weighting Function: no distance weighting
##
## Save grid: Salida/Tpi...
# Visualizando el mapa de índice de humedad topográfica
raster('Salida/Twi.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',palette = cpt(pal = 'idv_blues'))
# Visualizando el mapa de índice de posición topográfica
raster('Salida/Tpi.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'cont',palette = cpt(pal = 'neota_base_yellow'))
rsaga.get.usage(lib = 'ta_channels',module = 5)
## Search for SAGA command line program and modules...
## Done
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_channels
## library : Channels
## Usage: saga_cmd ta_channels 5 [-DEM <str>] [-DIRECTION <str>] [-CONNECTION <str>] [-ORDER <str>] [-BASIN <str>] [-SEGMENTS <str>] [-BASINS <str>] [-NODES <str>] [-THRESHOLD <num>]
## -DEM:<str> Elevation
## Grid (input)
## -DIRECTION:<str> Flow Direction
## Grid (optional output)
## -CONNECTION:<str> Flow Connectivity
## Grid (optional output)
## -ORDER:<str> Strahler Order
## Grid (optional output)
## -BASIN:<str> Drainage Basins
## Grid (optional output)
## -SEGMENTS:<str> Channels
## Shapes (output)
## -BASINS:<str> Drainage Basins
## Shapes (output)
## -NODES:<str> Junctions
## Shapes (optional output)
## -THRESHOLD:<num> Threshold
## Integer
## Minimum: 1
rsaga.geoprocessor(lib = 'ta_channels',module = 5,
param = list(DEM = 'Insumos/DEM_fill.sdat',
ORDER = 'Salida/Strahler_order',
SEGMENTS = 'Salida/red_drenaje'),
env = env)
##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_channels
## library : Channels
## tool : Channel Network and Drainage Basins
## author : O.Conrad (c) 2003
## processors : 4 [4]
##
## loading: DEM_fill
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Elevation: DEM_fill
## Flow Direction: <not set>
## Flow Connectivity: <not set>
## Strahler Order: Strahler Order
## Drainage Basins: <not set>
## Channels: Channels
## Drainage Basins: Drainage Basins
## Junctions: <not set>
## Threshold: 5
##
## Flow Direction
## Stream Order
## Junctions
## Drainage Basins
## Vectorising Grid Classes
##
##
## Parameters
##
##
## Grid system: 30; 558x 661y; 963390.8124x 8058312.3304y
## Grid: Drainage Basins
## Polygons: Polygons
## Class Selection: all classes
## Vectorised class as...: one single (multi-)polygon object
## Keep Vertices on Straight Lines: no
##
## class identification
## Create index: Drainage Basins
## edge detection
## edge collection
## Channels
## Save grid: Salida/Strahler_order...
## Save shapes: Salida/red_drenaje...
raster('Salida/Strahler_order.sdat') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_raster(style = 'fixed',breaks = seq(0,5,1),palette = '-viridis')
read_sf('Salida/red_drenaje.shp') %>%
tm_shape() +
tm_basemap('Esri.WorldTopoMap') +
tm_lines()